template<class T>
class Poorman 
{
	public:

	Poorman(void)
	{	
		boost::math::normal_distribution<> *d=new boost::math::normal_distribution<>(0,1); 
		_d=d;
		_conv=0;
	}
	~Poorman()
	{
		delete _d;
	}
	void operator()(mat X,double *Y,mat m, double s,int M)
	{
		
		T TdG(s,X,Y,1);
		mat G1=TdG.Sample(M);
		_m=(mean(G1(span((int)(0.1*M),M-1),span::all),0)).t();
		_S=-inv(Hessian(X,Y,_m,s));
		_S=cov(G1(span((int)(0.1*M),M-1),span::all));

	}

	inline double abserr(mat X)
	{
		double sum=0;
		int n=X.n_rows;
		for(int i=0;i<n;i++)
		{
			sum+=abs<double>(X(i,0));
		}
		return sum;
	}

	inline	mat Hessian(mat X,double *Y,mat theta,double s)
	{
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			
			xb=as_scalar(X(i,span::all)*theta);
			double pxb=pdf(*_d,xb);
			double gpxb=cdf(*_d,xb);
			mat xx=X(i,span::all).t()*X(i,span::all);	
			if(Y[i]==1)
			{	
				mat temp=(xb*pxb/gpxb+pow(pxb,2)/pow(gpxb,2))*xx;
				sum=sum+temp;				
			}else{
				mat temp=(-xb*pxb/(1-gpxb)+pow(pxb,2)/pow(1-gpxb,2))*xx;
				sum=sum+temp;				
			}
		}
		return -sum-1/s;

	}
	mat Get_Sig(void){return _S;}
	mat Get_Mu(void){return _m;}
	int Get_conv(void){return _conv;}
	private:
	mat _m;
	mat _S;
	int _conv;
	boost::math::normal_distribution<> *_d; 
};

